Thermophysical Properties and Phase Behavior of CO2 with Impurities: Insight from Molecular Simulations

Experimentally determining thermophysical properties for various compositions commonly found in CO2 transportation systems is extremely challenging. To overcome this challenge, we performed Monte Carlo (MC) and Molecular Dynamics (MD) simulations of CO2 rich mixtures to compute thermophysical properties such as densities, thermal expansion coefficients, isothermal compressibilities, heat capacities, Joule–Thomson coefficients, speed of sound, and viscosities at temperatures of (235–313) K and pressures of (20–200) bar. We computed thermophysical properties of pure CO2 and CO2 rich mixtures with N2, Ar, H2, and CH4 as impurities of (1–10) mol % and showed good agreement with available Equations of State (EoS). We showed that impurities decrease the values of thermal expansion coefficients, isothermal compressibilities, heat capacities, and Joule–Thomson coefficients in the gas phase, while these values increase in the liquid and supercritical phases. In contrast, impurities increase the value of speed of sound in the gas phase and decrease it in the liquid and supercritical phases. We present an extensive data set of thermophysical properties for CO2 rich mixtures with various impurities, which will help to design the safe and efficient operation of CO2 transportation systems.


INTRODUCTION
−3 The 2023 statistical review of world energy indicates that more than 80% of world energy consumption comes from fossil fuels. 2,4It is unfeasible to shift entirely to renewable energy resources on a very short time scale. 5Carbon Capture and Storage (CCS) is the most popular and large-scale process used in industries to meet the anthropogenic CO 2 emission targets. 6CCS is carried out in three different stages, namely, capture, transportation, and sequestration. 7Carbon capture and sequestration depends upon safe and economical transportation of CO 2 . 8In the past four decades, pipelines have been used to successfully inject pure CO 2 into depleted or nearly oil/gas fields for Enhanced Oil/gas Recovery (EOR). 3The injection of CO 2 captured from the flue gas stack is significantly different from pure CO 2 injection for EOR. 9 The difference is the presence of impurities, since no gas separation process is 100% efficient. 10mpure CO 2 also differs in composition depending on the source and technology of capture. 11It is possible to obtain pure CO 2 from impure CO 2 , but this will result in additional costs and carbon footprint. 11,12The most efficient and preferred phase of transporting CO 2 in the pipeline is the dense supercritical or liquid phase. 8The presence of impurities, especially noncondensable gases (for example, Ar, N 2 , H 2 , O 2 , and CH 4 which have a low boiling point compared to CO 2 ), reduces the density of impure CO 2 mixtures and is likely to introduce two-phase flow behavior. 13Two-phase flows during transportation lead to numerous challenges, mainly pressure surge, 14 which will sequentially lead to possible pipeline failure. 15The recommended level of major impurities in CO 2 rich stream for safe and efficient pipeline transportation and sequestration from different standards (NETL, 16 Dynamis, 17 and ISO 18 ) and projects (Porthos 19 and Teeside 13 ) are summarized in Tables 1 and 2. For quality standards of additional minor impurities in CO 2 transportation systems, the reader is referred to the review article of Simonsen et al. 18 CO 2 rich stream with impurities significantly alters the thermodynamic and transport properties of CO 2 , which will, in turn, impact the overall flow behavior, pipeline capacity, and operating window in CO 2 pipeline systems. 20,21Therefore, knowledge on thermodynamic and transport properties is indispensable to model the flow and phase behavior of impure CO 2 rich mixture within the operating window for the safe design and efficient operation of CO 2 transportation systems. 9Long-range tail corrections for Lennard-Jones (LJ) interactions are used for all components.
Table 2. CO 2 Quality Standards from the National Energy Technology Laboratory (NETL), 16 Dynamis, 17 International Standard Organization (ISO), 18 Porthos, 19 and Teesside The impurity percentages with an asterisk (*) in the NETL 16 advised limits imply that the total impurity concentration should be ≤4%.

Journal of Chemical & Engineering Data
impurity levels ranging from (1−10) mol % (which includes 12 CO 2 binary mixtures with 1 mol %, 5 mol %, and 10 mol % impurities which are shown in Section S14 of the Supporting Information and 36 multicomponent CO 2 mixtures with impurities ≤4 mol % which are shown in Sections S17 and S18 of the Supporting Information) using MC and MD simulation techniques.Simulations were carried out for a range of temperatures from −20 to 40 °C and pressures up to 200 bar.The main impurities, N 2 , Ar, H 2 , and CH 4 , are selected to investigate the effect of impurities on the thermodynamic and transport properties.A comprehensive list of chemical components, CAS numbers, and force fields of all components used in this work is shown in Table 1.Properties calculated within the operating window include density, isothermal compressibility, thermal expansion coefficient, heat capacity at constant volume and pressure, Joule−Thomson coefficient, shear viscosity, and speed of sound.We showed that, in comparison to pure CO 2 density, a CO 2 rich mixture with molecular weight lower than pure CO 2 at a condition had lower densities.We also showed that impurities decrease the value of thermal expansion coefficients, isothermal compressibilities, heat capacities, and Joule−Thomson coefficients in the gas phase and increase the value of these properties in the liquid and supercritical phases.Conversely, impurities tend to increase the speed of sound in the gas phase and decrease the speed of sound in the liquid and supercritical phases.Our results show that the order of influence of a particular impurity on a thermodynamic property other than density correlates with the critical temperature of that impurity.This article is organized as follows: The methodology used to compute thermodynamics and transport properties is explained in Section 2, followed by the simulation details in Section 3. In Section 4, we present the validation and results of the computed thermodynamic and transport properties.In Section 5, our key findings are summarized.

THEORETICAL BACKGROUND
To compute the speed of sound (c), one requires other properties, which include heat capacity at constant pressure (C P ), heat capacity at constant volume (C V ), and isothermal compressibility (β T ) 30,31

= c T P vC T P MC T V T P ( , )
( , ) ( , ) ( , ) where v is the molar volume of the pure component or the mixture and M is the molar mass of the pure component or the mixture.For a mixture, M can be calculated from the pure component molar mass where n is the number of components present in the mixture and x i and M i are the mole fraction and molar mass of each component present in the mixture.To calculate C V , C P , and β T , the derivatives of internal energy, volume, and enthalpy with respect to temperature and pressure have to be determined 36,45 i k j j j y where U and H are the internal energy and enthalpy of the system, respectively, and ⟨•••⟩ denotes the ensemble average of an ensemble.The internal energy (U = U internal + U external ) and enthalpy (H = U internal + U external + K + PV) in eq 3 and eq 4 include the kinetic energy term (K) in addition to the potential energy contribution from intramolecular molecular interaction (interaction inside molecules which is denoted as U internal ) and intermolecular interactions (interaction between molecules which is denoted as U external ).Hence, C V and C P have been split into ideal and residual contributions.Following the work of Lagache et al. 45 we can write i k j j j j j y { z z z z z i k j j j j j y { z z z z z Heat capacities be obtained from the standard thermodynamic databases 46,47 or from quantum mechanical calculations.In this study, the latter approach was used to calculate C V ideal (T) and C P ideal (T) using the Gaussian 09 software 48 where k B is the Boltzmann constant, N is the number of molecules in the system, U ̂(U ̂= U internal + U external ) is the configurational energy, and H ̂(H ̂= U internal + U internal + PV) is the configurational enthalpy of the system.The molar heat capacities (c P and c V ) can be obtained from where N A is Avogadro's number and N is the number of molecules in the system.The derivative ( ) is required to calculate β T defined in eq 5, and is computed using the following fluctuation formula, 49 i k j j j y Analogous to the computation of the speed of sound, the calculation of the Joule−Thompson coefficient requires heat capacity at constant pressure (c P ) and thermal expansion coefficient (α P ).
The value of α P is computed from the derivative of volume with respect to temperature i k j j j y The derivative ( ) which is required to calculate α T is computed using 36,42,45 It is important to note that C V and α P can be computed indirectly using the thermodynamic relations 50 The proofs of the mathematical derivation of equations to compute C P , C V , β T , c, α P , and μ JT from simulations are provided in Sections S3−S8 of the Supporting Information.

SIMULATION DETAILS
The Lorentz−Berthelot mixing rules were used to compute LJ interactions for dissimilar atoms: 49,57 Different force fields for CO 2 are available in the literature.These range from simple single-site force fields like the Higashi model 58 and Statistical Associating Fluid Theory (SAFT)-γ 59 to more complex three-site force fields such as the Transferable Potentials for Phase Equilibria (TraPPE)-rigid, 60 TraPPEflex, 61 Elementary Physical Model 2 (EPM2), 62 Zhang model, 63 and Cygan model. 64Aimoli et al. 42 65 Murthy et al. model, 66 and Galassi and Tildesley model. 67Force fields from Galassi and Tildesley 67 and Murthy et al. 66 were not optimized for VLE calculations. 60The force fields TraPPE 60 and Koster et al. 65 are three-site models with quadrupole moment, but they differ in their parametrization.Rahbari et al. 36 compared the performance of five different force fields for H 2 (Cracknell, 68 Buch, 69 Hirschfelder et al., 70 Marx and Nielaba, 71 and Koster et al. 65 ).These authors found that the force field developed by Koster et al. 65 best predicted the second derivative thermodynamic properties with the least deviation compared to data obtained from REFPROP. 72The TraPPE 60,73 force field was used in this study for CO 2 , N 2 , and CH 4 molecules since the primary objective in developing the TraPPE force field is to predict the thermophysical properties for a wide range of state conditions and compositions.Hydrogen is simulated using the Koster et al. 65 model.The force field from García-Pérez et al. 74 75−77 In GE, two simulation boxes are considered: one represents the liquid phase, and the other represents the gas phase.The simulation boxes are allowed to exchange energy, volume, and molecules.In a dense liquid simulation box, the insertion of a molecule in a single step is impeded due to the low probability of finding a cavity to accommodate a molecule, and the deletion of a molecule in a single step leaves the simulation box with a high energy penalty to form a new configuration. 53The CFCMC method 51−55 overcomes this drawback by gradual insertion and removal of so-called fractional molecules by which the surrounding whole molecules can adapt simultaneously by performing trial moves related to fractional molecules besides thermalization trial moves such as translations, rotations, and, volume changes.In CFCGE, two simulation boxes with indistinguishable whole molecules and fractional molecules (fractional molecules can be in either of the simulation boxes, but one per component type) are used for simulating the phase coexistence.The interaction of the distinguishable fractional molecule of a component type i with a whole molecule is scaled with a coupling parameter λ i ∈ [0, 1]. 78The trial moves related to fractional molecules are randomly changing the value of λ while keeping the orientation and position of all the molecules constant, insertion of a fractional molecule in another simulation box at a randomly selected orientation and position while keeping the orientation and position of the whole molecules constant, and changing the identity of a fractional molecule in a simulation box to a whole molecule while simultaneously transforming the randomly selected whole molecule to a fraction molecule in another simulation box. 79Further details specific to phase equilibria calculation in the CFCGE can be found elsewhere. 53,78The VLE of pure component (unary) systems, CO 2 , Ar, N 2 , H 2 , and CH 4 , were computed in the NVT version of the CFCGE.The phase equilibria (Pxy diagram) of binary mixtures which include CO 2 /Ar, CO 2 /N 2 , CO 2 /H 2 , and CO 2 /CH 4 are computed in the NPT version of the CFCGE.The simulation box sizes and initial distribution of whole molecules between the simulation boxes were specified based on the experimental data of the state point for both the unary and binary systems.For example, for phase equilibria computation of the CO 2 and Ar binary system at 105 bar, we chose 360 CO 2 molecules and 140 Ar molecules for the intended liquid box of size 30 Å and 180 CO 2 molecules and 320 Ar molecules for the intended gas box of size 40 Å.For a pressure of 35 bar, we chose 470 CO 2 molecules and 30 Ar molecules for the intended liquid box of size 30 Å and 325 CO 2 molecules and 175 Ar molecules for the intended gas box of size 40 Å.For simulating the phase coexistence of unary and binary systems, an equilibration run of 5 × 10 4 and 1 × 10 5 MC cycles was performed, respectively.The number of trial moves in an MC cycle in Brick-CFCMC equals the total number of molecules in the simulation box, with a minimum of 20.Following the equilibration run, a production run of 1 × 10 5 cycles was performed for unary systems and 2 × 10 5 cycles for binary systems to compute the coexistence densities and mole fractions of the components, respectively.
To calculate thermodynamic properties, density (ρ), isothermal compressibility (β T ), thermal expansion coefficient (α P ), molar heat capacity at constant volume (c V ), molar heat capacity at constant pressure (c P ), Joule−Thompson coefficient (μ JT ), and speed of sound (c) MC simulations were performed without fractional molecules.The values of c P , β T , μ JT , and α P were computed using the fluctuation equations in the NPT ensemble defined in eqs 4, 5, 12, and 13, respectively.To compute c V and consequently the speed of sound (c), it is important to extract the ensemble-averaged volume (⟨V⟩) from the NPT ensemble that reflects the same state for performing a simulation in the NVT ensemble with the same number of molecules.The simulations were performed with 300 molecules irrespective of the state and system.An equilibration of 5 × 10 4 MC cycles is performed to equilibrate the system successively, and 10 6 production runs are performed for each simulation.
In MD simulations, initial configurations of molecules in the cubic simulation boxes are constructed using PACKMOL 80 and fftool. 81Periodic boundary conditions are applied to simulation boxes in all directions.A cutoff radius of 12 Å is used for LJ interactions with analytic tail corrections.All molecules are treated as rigid bodies, and Newton's equations of motion are integrated using the velocity-Verlet algorithm with a time step of 0.5 fs.To thermostat and barostat the system, the Nose−Hoover type is used with coupling constants of 0.1 and 1 ps, respectively.The Particle−Particle Particle− Mesh (PPPM) method is used to handle long-range electrostatic interactions with a cutoff radius of 12 Å and 10 −6 accuracy.Shear viscosities (η) are calculated by performing an Equilibrium Molecular Dynamics (EMD) simulation by using the On-the-fly Computation of Transport Properties (OCTP) plugin in LAMMPS.The OCTP combines the Einstein relations with an order-n algorithm to calculate viscosity.Additional details about the OCTP's computation of transport properties can be found elsewhere. 82System sizes chosen to evaluate η of pure and multicomponent mixtures were 100 molecules for the very diluted gas phase, 300 molecules for the gas phase (40 to 80 bar), and 400 molecules for the liquid and supercritical phases (≥100 bar).The various steps involved in computing η for one independent simulation are as follows: First, a simulation is carried out in the NPT ensemble (0.5 ns equilibration run and 1 ns production run) to compute the ensemble average volume (⟨V⟩).Next, the simulation box is scaled according to the computed value of ⟨V⟩, and this system is then used to perform a simulation in the NVT ensemble (0.5 ns equilibration run and 1 ns production run) to compute the average total energy of the system (⟨E⟩).Finally, the ensemble average total energy is used to scale to the kinetic energy of the system to perform simulations in the NVE ensemble.Viscosities are calculated in the NVE ensemble, ensuring that the thermostat and barostat have no effect on the results.In the NVE ensemble, a production run of 5 ns is simulated to compute η.

RESULTS AND DISCUSSION
This section shows and discusses the thermodynamic and transport properties of single-component and multicomponent CO 2 systems.All properties of interest computed in this work from molecular simulations are compared with data sets generated from the NIST REFPROP database version-10.0 72 except when stated otherwise.EoS and correlation models used by REFPROP for computing pure component (CO 2 , N 2 , Ar, H 2 , and CH 4 ) thermodynamic and transport properties are listed in Table S11 of the Supporting Information.For multicomponent systems, the Groupe Europeèn de Recherches Gazieres (GERG-2008) EoS 24 is chosen in REFPROP.The GERG-2008 EoS is less accurate for unary systems. 72herefore, it is used only to compute the thermodynamic properties of multicomponent systems.The GERG-2008 EoS was originally developed for natural gas mixtures containing 21 components, including CO 2 , N 2 , Ar, H 2 , and CH 4 .Multicomponent mixtures containing high CO 2 concentrations and low levels of impurities were not the main focus for the development of the GERG-2008 EoS. 24,83The quality and quantity of experimental data used in developing the GERG-2008 EoS limit its accuracy. 24Thermophysical properties computed from molecular simulations will help in optimizing the EoS. 26Viscosities computed from MD simulations of multicomponent systems are compared with those obtained from the Extending Corresponding States (ECS) model available in REFPROP. 84Further details specific to the ECS model can be found elsewhere. 84Since the experimental data for multicomponent systems are scarce and simulating an exact composition as experiments is impossible with a small system of only 300 molecules, we opted to compare results from molecular simulations with those obtained from EoS.The numerical data used to generate all plots is provided in Sections S3, S4, and S14−S17 of the Supporting Information.The deviations of properties computed from simulations with respect to REFPROP data sets are computed using where χ Simulation  Critical temperatures (T c ) and densities (ρ c ) of all pure components computed from simulations using the law of rectilinear diameters and EoS models are listed in Table 3.Our results show that the computed liquid and vapor densities are in excellent agreement with respect to the REFPROP data set, except for the liquid densities of hydrogen.Deviations of the computed liquid densities and, in turn, the T c and ρ c of hydrogen with respect to REFPROP are due to the domination of quantum effects at low temperatures. 65This work focuses on temperatures significantly higher (>250 K) than H 2 VLE temperatures; hence, only the gas densities were considered for the validation of H 2 force field.
Figure 1 shows the phase equilibria (Pxy diagram) of four binary mixtures comprised of CO 2 , Ar, N 2 , H 2 , and CH 4 computed from CFCGE simulations (in the NPT version) compared with experimental data and the GERG-2008 EoS.The Pxy diagram, i.e., the bubble and dew points of CO 2 /Ar mixtures, is computed at 253.28 K and compared with experimental data of Coquelet et al., 85 and for CO 2 /CH 4 mixtures, simulations were performed at 250 K and compared to experimental data of Wei et al. 86 and Davalos et al. 87 Similarly, the Pxy diagrams of CO 2 /H 2 and CO 2 /N 2 mixtures are computed at 250 K and validated with experimental data of Tsang and Street 88 and Brown et al., 89 respectively.Figure 1a shows CO 2 liquid and gas mole fractions of CO 2 /Ar mixtures, computed from the CFCGE simulations compared to experimental data of Coquelet et al. 85 and the GERG-2008 EoS.Our results match well with the GERG-2008 EoS compared to the experimental data of Coquelet et al. 85 Bubble and dew points of CO 2 /CH 4 mixtures (shown in Figure 1b) computed from simulations agree well with experimental data of Davalos et al. 87 when compared to experimental data of Wei et al. 86 and the GERG-2008 EoS at high pressures.The phase equilibria of the CO 2 /H 2 binary system are shown in Figure 1c.REFPROP fails to converge for pressures larger than 170 bar when using the GERG-2008 EoS.This discrepancy seen in the GERG-2008 EoS was also reported in the work of Shin et al. 34 The bubble points obtained from the GERG-2008 EoS using REFPROP provide poor estimates when compared to the experimental data of Tsang and Street. 88Mole fractions computed from simulations agree fairly with experimental data of Tsang and Street, 88 and less than 5% relative deviation was noticed at high pressures.Dew points of CO 2 /N 2 mixtures shown in Figure 1d have a reasonable agreement with EoS and experimental data of Brown et al., 89 but the computed bubble points have a maximum relative deviation of 3% at high pressures.At high pressures, the dew and bubble points computed from CFCGE simulations agree moderately with the GERG-2008 EoS and experiments for all systems shown in Figure 1.This is because the mixing rules used in the simulations did not include binary interaction parameters.Computing bubble and dew points close to the critical pressure is challenging, even with a system size of 1000 molecules (the system size chosen in this study to compute the phase equilibria).In principle, one can perform simulations with large system sizes for pressures in the neighborhood of the critical pressure at the cost of much larger computations.Nevertheless, it is impossible to compute accurate bubble and dew points very close to critical pressure since simulation boxes may switch identities, which complicates ensemble averaging.Similarly, the GERG-2008 EoS model fails to predict reliable results for pressures close to the critical pressure due to the unavailability of experimental data to develop the GERG-2008 EoS.   The simulations are performed at 253.28 K for CO 2 /Ar mixtures and at 250 K for CO 2 /CH 4 , CO 2 /H 2 , and CO 2 /N 2 mixtures.

Thermal Expansion Coefficients (α P )
. Thermal expansion coefficients (α P ) computed from MC simulations using eq 13 in the NPT ensemble for pure CO 2 , binary CO 2 mixtures, and ternary CO 2 rich mixtures with different impurities are shown as a function of temperature and pressure in Figure 3. α P values of pure CO 2 computed at temperatures of 253, 273, 293, and 313 K are shown in Figure 3a.α P of pure CO 2 increases with temperature, whereas, with increasing pressure, the value of α P increases until it reaches a maximum value close to its saturation pressure and then decreases with increasing pressure.MC simulations closely predicted the peak of α P .Thermal expansion coefficients of pure CO 2 computed at conditions close to its saturation pressure were found to have large uncertainties and deviations when compared to the Span and Wagner EoS.Excluding conditions close to the critical point, α P of pure CO 2 computed from MC simulations was in good agreement with the Span and Wagner EoS for all temperatures with a maximum relative deviation of ca.6.8% at 293 K and 100 bar.In binary mixtures containing 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ), a similar pattern of α P with a peak value near its saturation pressure was found, as shown in Figure 3c at temperatures of 253 and 313 K.The computed α P of binary mixtures CO 2 /N 2 , CO 2 /CH 4 , CO 2 /Ar, and CO 2 /H 2 agrees fairly with α P obtained from the GERG-2008 EoS and has a maximum relative deviation of ca.13.9% at 313 K and 120 bar, ca.10.2% at 313 K and 60 bar, ca.8.1% at 313 K and 60 bar, and ca.20.4% at 313 K and 120 bar, respectively, excluding state points close to the critical point.State points of binary mixtures close to the critical point were identified based on large uncertainties observed in simulations and large relative deviation with the GERG-2008 EoS.The uncertainties and relative deviation of CO 2 /H 2 at 313 K decrease with increasing pressure.This suggests that the maximum relative deviation for the CO 2 /H 2 mixture observed in simulations at 313 K and 120 bar may be near its saturation/Widom line.The binary mixture containing 5 mol % H 2 increases the value of α P the most in liquid and supercritical phases when compared to α P of pure CO 2 , followed by binary mixtures containing 5 mol % N 2 , Ar, and CH 4 .In the gas phase, the presence of H 2 as an impurity decreases the value of α P followed by CH 4 , N 2 , and Ar compared to α P of pure CO 2 .Thermal expansion coefficients were also computed for 1 mol % and 10 mol % impurities at temperatures of 253, 273, 293, and 313 K from MC simulations.The computed thermal expansion coefficients provided in Tables S73−S167 of the Supporting Information show that the value of α P is affected based on impurities and concentration level of impurities.Thermal expansion coefficients of ternary mixtures in the liquid phase with 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ) are shown in Figure 3b.The effect of a particular mixture on α P was not discernible since the thermal expansion coefficients of all of the ternary mixtures computed at 253 K were within the limits of computed uncertainty of other ternary mixtures.Similarly, no significant differences were found between the thermal expansion coefficients of ternary mixtures obtained from the GERG-2008 EoS.

Isothermal Compressibilities (β T ).
The isothermal compressibility (β T ) computed from MC simulations using eq 5 in the NPT ensemble for pure CO 2 , binary CO 2 mixtures, and ternary CO 2 rich mixtures with different impurities are shown as a function of temperature and pressure in Figure 4. β T values of pure CO 2 at temperatures of 253, 293, 293, and 313 K are shown in Figure 4a.The maximum value of β T was observed in the gas phase.This implies that the volume change rate in response to the change in pressure is maximum when the fluid acts like an ideal gas.β T rapidly decreases with an increase in pressure for pressures lower than the saturation pressure of a particular temperature.For pressures away from the saturation/Widom line, the change in the value of β T was insignificant for all temperatures, as seen in Figure 4a.shown in Figure 4c.Isothermal compressibilities of binary mixtures computed from MC simulations agree fairly with the GERG-2008 EoS.The maximum relative deviations of CO 2 / N 2 , CO 2 /CH 4 , CO 2 /Ar, and CO 2 /H 2 binary mixtures were ca.14.1% at 253 and 80 bar, ca.17.5% at 313 and 120 bar, ca.12.9% at 253 and 120 bar, and ca.22.5% at 313 and 140 bar, respectively, excluding state points close to the critical point.
Figure 4c shows that the presence of impurity increases β T in the liquid and supercritical phases.In contrast, the presence of an impurity decreases β T in the gas phase.Comparing the β T of all binary mixtures shown in Figure 4c, a binary mixture with H 2 as an impurity increases the value of β T the most in the liquid and supercritical phases, followed by N 2 , Ar, and CH 4 .In the gas phase, a binary mixture with H 2 as an impurity decreases the value of β T followed by the values for N 2 , Ar, and CH 4 .This pattern of β T with respect to phases remained the same for both higher (10 mol %) and lower (1 mol %) concentrations of impurities, which is provided in Tables S73− S167 of the Supporting Information.The liquid phase β T of ternary mixtures with 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ) at 253 K are shown in Figure 4b.From Figure 4b, it is obvious that the presence of impurities tends to increase the value of β T .Similar to thermal expansion coefficients, the effect of a particular type of impurity on β T was unclear since the change in the value of β T due to the presence of impurities was limited.However, β T computed from MC simulations for CO 2 rich ternary mixtures with different impurities agrees fairly with the GERG-2008 EoS.For instance, the CO 2 rich ternary mixture with Ar and CH 4 as impurities results in the smallest increase in β T , while the CO 2 rich ternary mixture with N 2 and H 2 as impurities leads to the largest increase in β T .The difference in the value of β T due to the presence of different impurity combinations estimated from MC simulations was consistent with the GERG-2008 EoS.
4.1.5.Isobaric and Isochoric Heat Capacities.The constant pressure and volume heat capacities are computed using eqs 3 and 4, respectively.Residual heat capacities C P residual and C V residual are computed from MC simulations in the NPT and NVT ensemble, respectively, to obtain c P and c V .Figure 5 shows the c P of pure CO 2 and binary CO 2 mixtures and ternary CO 2 mixtures as a function of temperature and pressure.c P of pure CO 2 computed from MC simulations compared to the Span and Wagner EoS at temperatures of 253, 293, 293, and 313 K are shown in Figure 5a.Heat capacities at a constant pressure of pure CO 2 are sensitive to pressures closer to the saturation/Widom line, where a sudden surge in c P is noticed.Figure 5a also shows that c P increases with temperature.The maximum value of c P observed in Figure 5a at 313 K is calculated by using the Span and Wagner EoS.However, this peak value of c P cannot be predicted by MC simulations unless one considers an extremely large system.This is because the correlation length, which measures the spatial extent of spontaneous density fluctuations, diverges near the critical point. 57c P of pure CO 2 computed from MC simulations resulted in minor systematic overprediction with a maximum relative deviation of ca.6.9% at 313 K and 120 bar.At 20 bar, relative deviations of c P with the Span and Wagner EoS are less than 1% for all temperatures shown in Figure 5a, except at 253 K.At 253 K and 20 bar, c P has a relative deviation of ca.3.8%, which is attributed to the close distance to the vapor−liquid phase transition.c P of binary mixtures with 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ) at temperatures of 253 and 313 K are shown in Figure 5b.Similar to c P of pure CO 2 , systematic overprediction of c P is observed in binary mixtures when compared to the GERG -2008 EoS.The maximum relative deviations of CO 2 /N 2 , CO 2 /CH 4 , CO 2 /Ar, and CO 2 /H 2 binary mixtures were ca.8.4% at 313 and 140 bar, ca.7.3% at 313 and 160 bar, ca.6.9% at 313 and 180 bar, and ca.15.7% at 313 and 140 bar, respectively, excluding state points close to the saturation/Widom line.Although c P values of different binary mixtures may appear to be identical in Figure 5b, impurities altered the value of heat capacities significantly compared to c P of pure CO 2 .The presence of impurities tends to increase the value of c P for pressures larger than the pressure at which the heat capacity peaks.Conversely, when the pressure is lower than the pressure at which the heat capacity peaks, the presence of impurities tends to decrease the value of c P .For instance, c P values of the CO 2 /H 2 binary mixture at 313 K and 160 bar were ca.13.2% larger than c P of pure CO 2 .c P of the CO 2 /H 2 binary mixture at 313 K and 60 bar were ca.13.7% smaller than c P of pure CO 2 .At conditions away from the saturation/ Widom line, a particular type of impurity did not influence c P to a great extent.However, at conditions close to the saturation/Widom line, the impact of impurities on c P becomes significant based on the type of impurity where a binary mixture with H 2 as an impurity increases c P to a great extent followed by N 2 , Ar, and CH 4 .Similar to binary mixtures with 5 mol % impurities, binary mixtures with 1 and 10 mol % impurities decreased and increased the value of c P compared to c P of pure CO 2 before and after the maximum heat capacity, respectively; see Tables S73−S167 of the Supporting Information.c P of ternary mixtures with 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ) at 253 K are shown in Figure 5c.From Figure 5c, it is clear that the presence of impurities tends to decrease the c P at 253 K. Figure 5c also confirms that the effect of a particular type of impurity combination on c P is insignificant.
Heat capacities at constant volume of pure CO 2 at temperatures of 253, 293, 293, and 313 K are compared with the Span and Wagner EoS in Figure 6a.The computed c V is in excellent agreement with c V obtained from the Span and Wagner EoS with a maximum relative deviation of ca. 5% at 313 K and 120 bar considering conditions away from the saturation/Widom line.c V values of pure CO 2 vary insignificantly for pressures larger than its saturation pressure irrespective of the temperature.Similarly, c V of binary mixtures with 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ) at temperatures of 253 and 313 K seen in Figure 6c have no significant difference in the value of c V between 253 and 313 K for pressures larger than its saturation pressure.The presence of an impurity tends to decrease c V hardly, irrespective of the type of impurity.c V of ternary mixtures with 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ) at 253 K are shown in Figure 6b.In the reduced axis range seen in Figure 6c, it is clear that the presence of impurities (in CO 2 /Ar/H 2 ) decreases the c V maximum by ca.3%.It is important to mention that c V can be calculated from c P using eq 15.However, the indirect computation of c V using eq 15 is prone to high statistical uncertainties, 90,91 which will subsequently affect the calculations of speed of sound.Hence, constant volume heat capacities are computed by sampling eq 3 in the NVT ensemble with a volume of the state obtained from the NPT simulation.
4.1.6.Joule−Thomson Coefficients (μ JT ).Joule−Thomson coefficients (μ JT ) of pure CO 2 , binary CO 2 mixtures, and ternary CO 2 rich mixtures computed using eq 12 as a function of the temperature and pressure are shown in Figure 7. μ JT decreases as the temperature increases in the gas phase.Conversely, μ JT increases with the temperature in the liquid and supercritical phases.The computed μ JT of pure CO 2 at temperatures of 253, 293, 293, and 313 K seen in Figure 7a agrees decently with μ JT obtained from the Span and Wagner EoS.The deviations of μ JT in the gas phase are larger when compared to the liquid and supercritical phases, and uncertainties of μ JT are quite significant for extremely low pressures.This is because of eq 12 used for computing μ JT .Joule−Thomson coefficients are computed indirectly using eq 12, and its uncertainties are computed using eq S84 presented in Section S11 of the Supporting Information.The uncertainties computed using eq S84 depend on the precision and accuracy of α P and c P .In particular, computation of μ JT is highly sensitive to α P .Even an insignificant error in the α P computation might result in high uncertainties and deviations in μ JT .A similar drawback of using eq 12 is also reported in ref 92.μ JT of binary mixtures with 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ) at temperatures of 253 and 313 K are shown in Figure 7b.The value of μ JT decreases due to the presence of impurities in the gas phase and increases in the liquid and supercritical phases.H 2 has the most significant effect on the values of μ JT , followed by those of N 2 , Ar, and CH 4 .μ JT computed from MC simulations for binary mixtures rich in CO 2 with 1 mol %, 5 mol %, and 10 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) at 253, 273, 293, and 313 K listed in Tables S73−S167 of the Supporting Information also shows the same qualitative trend with respect to impurity type.μ JT listed in Tables S73−S167 of the Supporting Information also indicates that the impact of impurities increases with increasing concentration of impurities.μ JT of ternary mixtures with 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ) at 253 K are shown in Figure 7c.μ JT of ternary mixtures computed from MC simulations at 253 K agrees well with the μ JT obtained from the GERG-2008 EoS, considering the uncertainty range.Figure 7c also shows that the presence of impurities tends to increase the μ JT of ternary mixtures significantly.The impact of impurities is crucial and should be considered in the initial computations, especially at conditions where the inversion of μ JT takes place.For instance, μ JT inversion occurs close to 80 bar for pure CO 2 at 253 K, but for impure CO 2 ternary mixtures at 253 K, μ JT inversion occurs at pressures larger than 120 bar.
4.1.7.Speed of Sound (c).The speed of sound (c) of pure CO 2 , binary CO 2 mixtures, and ternary CO 2 rich mixtures computed as a function of temperature and pressure using eq 1 from MC simulations is shown in Figure 8.Since c is a function of ρ, β T , c P , and c V , uncertainties associated with c are computed using eq S78 derived in Section S11 of the Supporting Information.The speed of sound of pure CO 2 computed at temperatures of 253, 293, 293, and 313 K results in overprediction compared to c obtained using the Span and Wagner EoS.Since c depends on ρ, β T , c P , and c V , small deviations in these dependent properties will overestimate c.The overestimation of c is majorly due to the minor overestimation of c P and underestimation of β T .For example, c of pure CO 2 is overestimated by a relative deviation of ca.6.7% at 253 K and 200 bar, while the overestimation of c P is ca.3.6% and the underestimation of β T is ca.8.4%.The magnitude of the overprediction was found to be less than 10% for all temperatures and pressures.Figure 8a shows that c decreases marginally until the pressure reaches its critical pressure and increases significantly for the pressure larger than the saturation pressure for a particular temperature.Figure 8a also shows that c is temperature dependent and decreases with increasing temperature.The impacts of impurities on c are evaluated by comparing c computed from MC simulation and the GERG-2008 EoS for binary mixtures consisting of 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ) as shown in Figure 8b.In the liquid and supercritical phases, the presence of H 2 as an impurity is observed to have the largest impact on decreasing the value of c, followed by those of CH 4 , Ar, and N 2 .In the gas phase, the presence of H 2 as an impurity is observed to have the largest impact on increasing the value of c followed by that of CH 4 , Ar, and N 2 .The computed c of binary mixtures rich in CO 2 with 1 mol % and 10 mol % impurity (CH 4 , Ar, N 2 , and H 2 ) for temperatures of 253, 273, 293, and 313 K from MC simulations and the GERG-2008 EoS, as listed in Tables S73−S167 of the Supporting Information, shows the same behavior with respect to type of impurity.Speeds of sound of ternary mixtures with 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ) at 253 K are shown in Figure 8c.MC simulations overpredicted the value of c compared to the GERG-2008 EoS.However, MC simulations closely predicted the decrease in the value of c due to the presence of impurities.For instance, the difference in the value of c at 200 bar between pure CO 2 and a ternary CO 2 mixture with N 2 and H 2 as impurities is ca.56 m/s from MC simulations and ca.46 m/s from the GERG-2008 EoS.
4.2.Transport Property.4.2.1.Viscosities (η).Viscosities (η) of pure CO 2 , binary CO 2 mixtures, and ternary CO 2 rich mixtures computed from MD simulations as a function of temperature and pressure are shown in Figure 9. Viscosities of pure CO 2 were computed at temperatures of 253, 273, 293, and 313 K.The computed viscosities are compared with those calculated from the correlation of Laesecke et al. 25 Our results show a good agreement with the model for all temperatures as seen in Figure 9a with a maximum relative deviation of ca.14.6% in the gas phase (at 313 K and 40 bar), excluding state points close to the saturation/Widom line.η at pressures below the saturation pressure for a temperature remains relatively constant, but for pressures far away from the saturation/ Widom line, η increases with an increasing pressure at a constant temperature.η is observed to be highly dependent on temperature compared to pressure as seen in Figure 9a, where η decreases with increasing temperature.To analyze the effect of impurities, η of binary mixtures with 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MD simulations at 253 and 313 K are compared with viscosities obtained from REFPROP in Figure 9b.Computed η of binary CO 2 mixtures qualitatively correlate well with η obtained from REFPROP.The existence of an impurity in a CO 2 mixture tends to reduce η compared to η of pure CO 2 .Based on η computed from MD simulations and REFPROP, it is clear that a binary mixture with H 2 as an impurity decreases η the most in liquid and supercritical phases.The uncertainty range of computed viscosities from MD simulations made it difficult to interpret a particular type of impurity that impacts η the most next to H 2 .However, viscosities obtained from REFPROP indicate that a binary mixture with N 2 as an impurity decreases η after H 2 followed by CH 4 and Ar for pressure away from the saturation/Widom line at 253 and 313 K.In addition to viscosity data of binary mixtures with 5 mol % impurities, binary mixtures with 1 mol % and 10 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) at temperatures of 253, 273, 293, and 313 K computed from MD simulations are listed along with data obtained from REFPROP in Tables S73−S167 of the Supporting Information.To analyze the effect on η due to the presence of a particular combination of impurities in CO 2 rich mixtures, η of ternary mixtures were computed from MD simulations and compared with η obtained from REFPROP at 253 K, as shown in Figure 9c.η of all ternary mixtures shown in Figure 9c has a concentration of 96 mol % CO 2 and 2 mol % for each of the two impurities (CH 4 , Ar, N 2 , and H 2 ).Comparing η of ternary mixtures obtained from REFPROP, it is clear that mixtures with H 2 as an impurity reduce the liquid viscosities of CO 2 mixtures to a larger extent compared to mixtures without H 2 .A similar qualitative trend of η was also observed from simulations at 40, 60, and 80 bar.For higher pressures, considering the computed uncertainties, the effect of a particular combination of impurities was difficult to evaluate, with a marginal decrease in the liquid viscosities.Uncertainties associated with viscosities can be reduced by performing multiple independent simulations, but this will increase computational costs. 93We refrained from conducting additional simulations due to the marginal difference in liquid viscosities observed between ternary mixtures in Figure 9c.
In summary, we performed molecule simulations of CO 2 rich mixtures with N 2 , Ar, CH 4 , and H 2 as impurities using the force fields mentioned in Section S1 of the Supporting Information.The thermodynamic and transport properties were computed at temperatures of 253, 273, 239, and 313 K and pressures up to 200 bar for pure CO 2 and binary mixtures rich in CO 2 with 1 mol %, 5 mol %, and 10 mol % impurities.The computed thermodynamic and transport properties were found to be in good agreement with the EoS for pure and binary systems, except at conditions close to the saturation/ Widom line.The thermodynamic and transport properties were also computed for 24 ternary and 12 quaternary CO 2 rich mixtures for various concentrations of impurities listed in Tables S169 and S362 of the Supporting Information, respectively, for temperatures of 253, 273, 239, and 313 K and pressures up to 200 bar.Results of the thermodynamic and transport properties of ternary and quaternary mixtures are provided in Sections S16 and S17 of the Supporting Information, respectively.The thermodynamic and transport properties of multicomponent CO 2 mixtures were compared with the GERG-2008 EoS 24 and the ECS model, 84 respectively, and showed good agreement.We show that molecular simulations are a powerful tool to compute the thermodynamic and transport properties of multicomponent mixtures seen in the transportation of CO 2 with a smaller system of 300 molecules.These thermophysical properties will help in modeling and designing pipelines for CO 2 transportation, which will be the focus of further work.

CONCLUSIONS
In this study, the effect of impurities in CO 2 rich mixtures on the value of thermodynamic and transport properties such as densities, thermal expansion coefficients, isothermal compressibilities, heat capacities at constant pressure, heat capacities at constant volume, Joule−Thomson coefficients, speed of sound, and viscosities were investigated using molecular simulations.The CFCMC method was used to compute the VLE of pure components, such as CO 2 , CH 4 , Ar, N 2 , and H 2 , to validate force fields used in molecular simulations.−97 The phase equilibria of CO 2 /Ar, CO 2 /CH 4 , CO 2 / N 2 , and CO 2 /H 2 binary mixtures were also computed using the CFCMC method and compared with the GERG-2008 EoS 24 and data from the literature, 85−89 showing a good agreement.The thermodynamic and transport properties were computed for pure CO 2 and binary and ternary mixtures rich in CO 2 at temperatures of 253, 273, 293, and 313 K and for pressures ranging from 20 to 200 bar using MC and MD simulations.The computed thermodynamic and transport properties of pure CO 2 are in excellent agreement with corresponding values obtained from the Span and Wagner EoS. 23Thermodynamic and transport properties of CO 2 rich binary mixtures with 1 mol %, 5 mol %, and 10 mol % concentrations of noncondensable impurities such as CH 4 , Ar, N 2 , and H 2 computed from simulations showed a good agreement with the GERG-2008 EoS. 24The computed thermodynamic and transport properties of CO 2 rich ternary mixtures with various impurities were compared with the GERG-2008 EoS, 24 showing a good agreement.The effect of different types of impurities on specific thermodynamic and transport properties was evaluated.Our findings show that CO 2 rich mixtures with impurities have low densities compared to densities of pure CO 2 .The magnitude of reduction in densities of a CO 2 rich mixture depends strongly on the molecular weight of impurities present in a mixture.Mixtures with molecular weight lower than pure CO 2 were observed to have reduced densities compared to pure CO 2 .CO 2 rich mixtures containing H 2 as an impurity led to the most significant decrease in the value of thermal expansion coefficients, isothermal compressibilities, heat capacities at constant pressure, and Joule− Thomson coefficients followed by N 2 , Ar, and CH 4 in the gas phase.In the liquid and supercritical phases, the presence of H 2 as an impurity led to the most significant increase in the values of thermal expansion coefficients, isothermal compressibilities, heat capacities at constant pressure, and Joule− Thomson coefficients followed by N 2 , Ar, and CH 4 .In contrast, the presence of H 2 as an impurity in the CO 2 rich mixture increased the value of speed of sound in the gas phase and decreased it in the liquid and supercritical phases.The order of effect due to a particular impurity on thermal expansion coefficients, isothermal compressibilities, heat capacities at constant pressure, Joule−Thomson coefficients, and speed of sound correlates with the critical temperature of impurities.In our investigation of heat capacities at constant Journal of Chemical & Engineering Data volume, we found that the presence of impurities did not have a significant impact.Finally, differences in the values of viscosities in CO 2 rich mixtures due to the presence of impurities were evaluated.Our findings showed that mixtures containing H 2 as an impurity significantly reduced viscosities in liquid and supercritical phases.

■ ASSOCIATED CONTENT
with the B3LYP theory and a 6-31G(d,p) basis set.The derivatives U T residual and H T residual required to calculate C V and C P are computed from fluctuations in NVT and NPT ensembles, respectively30,31,36,45 83 4.1.2.Densities (ρ).Densities of pure CO 2 computed from MC and MD simulations for temperatures of 253, 273, 293, and 313 K and pressures ranging from 20 to 200 bar are shown in Figure 2a.MC and MD simulations have an excellent agreement with the Span and Wagner EoS for all temperatures.The computed densities from MC simulations have a maximum relative deviation of ca.0.93% at 313 K and 200 bar, and MD simulations have a maximum relative deviation of

Figure 2 .
Figure 2. Computed densities as a function of temperature and pressure.(a) The calculated densities of pure CO 2 from MC simulations (closed symbols), MD simulations (open symbols), and the Span and Wagner EoS 23 (solid lines) for temperatures of 253, 273, 293, and 313 K. (b) Densities of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with densities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures of 253 and 313 K. (c) Densities of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with densities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.

Figure 3 .
Figure 3. Computed thermal expansion coefficients as a function of temperature and pressure.(a) The calculated thermal expansion coefficients of pure CO 2 from MC simulations (closed symbols) and the Span and Wagner EoS 23 (solid lines) for temperatures of 253, 273, 293, and 313 K. (b) Thermal expansion coefficients of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with thermal expansion coefficients of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures of 253 and 313 K. (c) Thermal expansion coefficients of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with thermal expansion coefficients of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.
Figure 4a also shows that β T is dependent on temperature and β T increases with increasing temperature.β T of pure CO 2 computed from MC simulations agrees qualitatively with the Span and Wagner EoS with a maximum relative deviation of ca.19.9% at 293 K and 100 bar.β T values of binary mixtures with 95 mol % CO 2 and 5 mol % of one of the impurities (CH 4 , Ar, N 2 , and H 2 ) at temperatures of 253 and 313 K are

Figure 4 .
Figure 4. Computed isothermal compressibilities as a function of temperature and pressure.(a) The calculated isothermal compressibilities of pure CO 2 from MC simulations (closed symbols) and the Span and Wagner EoS 23 (solid lines) for temperatures of 253, 273, 293, and 313 K. (b) Isothermal compressibilities of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with isothermal compressibilities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures of 253 and 313 K. (c) Isothermal compressibilities of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with isothermal compressibilities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.

Figure 5 .
Figure 5. Computed isobaric heat capacities as a function of temperature and pressure.(a) shows the calculated isobaric heat capacities of pure CO 2 from MC simulations (closed symbols) and the Span and Wagner EoS 23 (solid lines) for temperatures: 253, 273, 293, and 313 K. (b) shows isobaric heat capacities of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with isobaric heat capacities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures: 253 and 313 K. (c) shows isobaric heat capacities of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with isobaric heat capacities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.

Figure 6 .
Figure 6.Computed isochoric heat capacities as a function of temperature and pressure.(a) The calculated isochoric heat capacities of pure CO 2 from MC simulations (closed symbols) and the Span and Wagner EoS 23 (solid lines) for temperatures of 253, 273, 293, and 313 K. (c) Isochoric heat capacities of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with isochoric heat capacities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures of 253 and 313 K. (b) Isochoric heat capacities of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with isochoric heat capacities of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.

Figure 7 .
Figure 7. Computed Joule−Thomson coefficients as a function of temperature and pressure.(a) The calculated Joule−Thomson coefficients of pure CO 2 from MC simulations (closed symbols) and the Span and Wagner EoS 23 (solid lines) for temperatures of 253, 273, 293, and 313 K. (b) Joule−Thomson coefficients of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with Joule−Thomson coefficients of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures of 253 and 313 K. (c) Joule−Thomson coefficients of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with Joule−Thomson coefficients of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.

Figure 8 .
Figure 8. Computed speed of sound as a function of temperature and pressure.(a) The calculated speed of sound of pure CO 2 from MC simulations (closed symbols) and the Span and Wagner EoS 23 (solid lines) for temperatures of 253, 273, 293, and 313 K. (b) The speed of sound of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with the speed of sound of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 for temperatures of 253 and 313 K. (c) The speed of sound of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MC simulations (closed symbols) and the GERG-2008 EoS 24 (lines) compared with the speed of sound of pure CO 2 computed from MC simulations and the Span and Wagner EoS 23 at 253 K.

Figure 9 .
Figure 9. Computed viscosities as a function of temperature and pressure.(a) The calculated viscosities of pure CO 2 computed from MD simulations (closed symbols) compared with the correlation of Laesecke et al. (2017) 25 (solid lines) for temperatures of 253, 273, 293, and 313 K. (b) Viscosities of binary mixtures with 95 mol % CO 2 and 5 mol % impurities (CH 4 , Ar, N 2 , and H 2 ) computed from MD simulations (closed symbols) and REFPROP72 (lines) compared with viscosities of pure CO 2 computed from MD simulations and REFPROP72 for temperatures of 253 and 313 K. (c) Viscosities of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of two components (CH 4 , Ar, N 2 , and H 2 ) computed from MD simulations (closed symbols) and REFPROP72 (lines) compared with viscosities of pure CO 2 computed from MD simulations and REFPROP 72 at 253 K.

Table 1 .
Description of All Components Used in This Work a 13 a

Journal of Chemical & Engineering Data calculations.
The VLE of pure components (CO 2 , Ar, N 2 , H 2 , and CH 4 ) computed from the CFCMC simulations in the NVT version are compared to the REFPROP database.The computed VLE curves compared to the REFPROP database are shown in Section S2 of the Supporting Information.
and χ REFPROP are the properties of interest computed from simulations and REFPROP, respectively.For the sake of clarity, plots for temperatures of 253 and 313 K are shown for binary systems, and plots for a temperature of 253 K are shown for multicomponent systems despite that properties of interest were estimated for four different temperatures: 253, 273, 293, and 313 K.4.1.Thermodynamic Properties.4.1.1.Phase Equilibria.

Table 3 .
Comparison of Critical Temperatures (T c at temperatures of 253, 273, 293, and 313 K from MC simulations and MD simulations.The computed densities are provided in Tables S73−S167 of the Supporting Information.Densities of all binary mixtures are shown in Section S15 of the Supporting Information; densities decrease with increasing mol % of impurity compared to densities of pure CO 2 irrespective of the temperature.Figure2cshows the liquid densities of ternary mixtures with 96 mol % CO 2 and 2 mol % impurities for each of the two components (CH 4 , Ar, N 2 , and H 2 ) at 253 K.The liquid densities of ternary mixtures may appear to have some deviation when compared to the GERG-2008 EoS due to the reduced axis range seen in Figure2c.The maximum relative deviation of ternary mixture liquid densities in Figure2cwas ca.0.78% for the CO 2 /Ar/N 2 and CO 2 /N 2 /CH 4 mixtures at 200 bar.The liquid densities of ternary mixtures were observed to have good agreement with A comprehensive analysis of the gas phase densities of binary mixtures with 5 mol % of one of the impurities found that impurities do not significantly change the gas phase densities.The influence of impurities on densities is consistent with the molecular weight of the mixtures.A mixture with H 2 as an impurity decreases the density to a large extent, followed by Ar, N 2 , and CH 4 .Densities were also computed for 1 and 10 mol % impurities the GERG-2008 EoS.Similar to binary mixtures, ternary mixtures with the least molecular weight tend to have lesser densities when compared to other ternary mixtures.For example, a ternary mixture with 2 mol % H 2 and 2 mol % CH 4 as impurities, which have the lowest molecular weight among other ternary mixtures, have the lowest densities when compared to densities of other ternary mixtures.